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1.  Overview 


Mesotek+  is  a  parallel  C++  mesoscale  polymer  modeling  program  that  uses  the  field-theoretic 
formalism  developed  by  Glenn  Fredrickson  and  coworkers  at  the  University  of  California,  Santa 
Barbara.  This  code  is  generalized  for  a  variety  of  non-closed-loop  polymer  architectures  such  as 
linear,  star,  and  branched,  and  it  can  be  applied  to  model  the  equilibrium  morphology  of  both 
neat  polymeric  systems  and  nanocomposites.  Two  different  field  theories,  self-consistent  field 
theory  and  complex  Langevin  theory,  are  implemented.  Self-consistent  field  theory  is 
appropriate  when  thennal  fluctuations  can  be  ignored,  such  as  in  dense,  high-molecular  weight 
polymer  melts.  Complex  Langevin  theory  should  be  used  for  polymer  solutions  in  the  dilute  or 
semidilute  regimes,  microemulsion  and  micellar  phases,  block  copolymers  near  the  order- 
disorder  transition,  and  blends  near  a  critical  point,  second-order  phase  transition,  or  weakly 
first-order  phase  transitions  (/).  Although  complex  Langevin  theory  can  also  be  used  to  look  at 
dense  polymer  melts,  it  is  much  more  computationally  expensive  than  self-consistent  field  theory 
and  should  be  used  only  when  necessary. 

When  the  neat  polymeric  systems  are  modeled,  two  routes  can  be  taken  to  reach  equilibrium: 
diffusive  dynamics  or  optimization.  Although  these  two  methods  transverse  different  pathways, 
they  should  both  reach  the  same  equilibrium  state.  Diffusive  dynamics  can  be  used  if  the  kinetic 
pathway  toward  equilibrium  is  important;  however,  the  results  should  be  analyzed  with  caution 
since  it  does  not  include  hydrodynamic  and  elastic  stress  couplings,  which  are  often  important 
(1 ).  The  optimization  pathway  occurs  through  fictitious  dynamics  that  are  not  physically 
realistic.  If  only  the  equilibrium  morphology  is  desired,  optimization  should  be  used  because  it 
allows  for  a  quicker  route  to  equilibrium  than  diffusive  dynamics. 

Nanocomposites  can  only  be  modeled  using  the  optimization  route  with  self-consistent  field 
theory.  In  this  method,  the  polymers  are  modeled  using  self-consistent  theory  for  fixed 
nanoparticles,  while  the  nanoparticles  are  moved  using  sampling  methods  such  as  Monte  Carlo 
or  Langevin  (2). 

This  report  presents  an  introduction  on  how  to  use  Mesotek+  (it  is  assumed  that  the  reader  has  a 
basic  knowledge  of  field-theoretic  methods).  If  the  reader  is  interested  in  the  underlying  theory 
used  in  this  software,  a  variety  of  articles  are  available  (1,  3). 
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2.  Defining  Your  System 


To  run  a  simulation,  it  is  necessary  to  produce  an  input  file  that  includes  your  system’s 
parameters.  The  minimum  information  needed  is  the  simulation  box  size  (discretization  and 
dimensions),  details  regarding  your  species  (type  and  interaction),  and  components  (architecture 
and  fractions). 

2.1  Box  Size 

To  parameterize  your  simulation  box,  it  is  necessary  to  specify  the  box  angles,  the  resolution  of 
each  dimension,  and  the  length  of  the  simulation  box  (see  figure  1).  This  group  is  rather 
straightforward.  Here,  the  rank  is  the  dimension  of  your  system.  In  a  one-dimensional 
simulation,  the  dimensions  not  specified  are  assumed  to  be  homogeneous  such  that  only  the 
lamellar  or  disordered  phases  are  observed.  For  a  two-dimensional  (2-D)  simulation,  the  third 
direction  that  is  not  explicitly  simulated  is  assumed  to  be  unifonn  such  that  three-dimensional 
(3-D)  spheres  will  appear  as  cylinders  in  2-D.  The  simulation  box  is  assumed  be  periodic  unless 
a  confined  geometry  is  specified.  The  parameters  n_i  and  d_i  are  the  number  of  collocation 
points  and  lattice  length,  respectively,  in  the  i  direction  (i  =  1,  2,  or  3).  If  complex  Langevin  is 
used,  it  is  also  necessary  to  include  the  line  fftwmode  =  c2c,  which  specifies  the  use  of  complex 


fields. 


2-direction 


n  1 


Note:  Box  size  (rank  =  3;  n_l  =  32;  n_2  =  64;  n_3  =  16;  d_l  =  8.0;  d_2  =  12.0;  d_3  =  4.0; 
ang_l_2  =  120;  ang_l_3  =  90;  ang_2_3  =  90). 


Figure  1.  Simulation  box  that  includes  the  parameters  defined  in  the  box  size  group. 
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2.2  Species 


The  species  indicate  the  species  in  your  system  as  well  as  the  interactions.  Mesotek+  has  two 
different  types  of  inputs:  experimental  and  theoretical. 

2.2.1  Theoretical 

When  theoretical  units  are  chosen,  all  length  scales  are  scaled  by  the  polymer’s  radius  of 
gyration  (Rg),  and  the  statistical  segment  length  is  assumed  to  be  equal  for  each  species.  A 
subgroup  must  be  produced  for  each  species  as  well  as  for  all  interactions.  For  example,  for  an 
AB  homopolymer  blend  or  any  other  polymer  architecture  with  two  polymer  species  (e.g.,  AB, 
ABA,  etc.),  you  must  make  a  subgroup  for  A  and  a  subgroup  for  B,  as  well  as  a  sub-subgroup  for 
the  interaction  chi  AB.  If  diffusive  dynamics  are  used  instead  of  optimization,  the  diffusion 
coefficient  must  also  be  specified  through  the  parameter  diffusivity.  The  sample  code  for  this 
model  is  in  code  la. 


Code  la:  Sample  Code  for  Theoretical 

Code  lb:  Sample  Code  for  Experimental 

Parameters 

Species) 

Species) 

A) 

A{ 

Rg2/M  =1.; 

Rg2/M  =1.; 

mass  density=  1.; 

mass  density=  1.; 

type  =  Gaussian; 

type  =  Gaussian; 

} 

diffusivity  =  1 ; 

B) 

} 

Rg2/M  =1.; 

mass  density=  1.; 

Rg2/M  =1.; 

type  =  Gaussian; 

mass  density=  1.; 

} 

type  =  Gaussian; 

C) 

diffusivity  =  1 ; 

Rg2/M  =1.; 

} 

mass  density  =1.; 

chi  { 

type  =  Solvent; 

A  B{ 

} 

~  A  =  0,2; 

chi  ) 

} 

A  B) 

} 

A  =  0.2; 

} 

} 

B  C) 

A  =  0.1; 

} 

A  C) 

A  =  0.2; 

} 

} 

} 
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If  an  additional  species  is  added,  such  as  a  solvent,  an  additional  species  group  C  would  have  to 
be  defined  along  with  two  additional  sub-subgroups  for  the  interactions  chi  AC  and  chi  BC  (see 
code  lb).  In  this  example,  the  interaction  between  the  different  species  is  given  by  the  variable 
A,  which  is  the  minimum  constant  required  to  define  X,  the  Flory-Huggins  interaction  parameter. 

2.2.2  Experimental 

When  experimental  units  are  chosen,  the  length  scales  are  in  physical  dimensions  (angstroms). 
The  code  structure  for  species  is  similar  to  the  structure  used  in  the  theoretical  system  but  will 
include  additional  information  about  the  experimental  system.  When  the  polymers  are  described, 
it  is  necessary  to  input  two  out  of  three  values  per  species  (mass  density,  molecular  weight  of  a 
bead  [MO],  or  molar  volume  of  a  bead)  as  well  as  one  experimental  length  scale  (either  the  end- 
to-end  distance  squared  per  molecular  weight  [R0A2/M]  or  the  radius  of  gyration  squared  per 
molecular  weight  [Rg2/M]).  The  end-to-end  distance  is  related  to  the  radius  of  gyration  by  a 
simple  constant  (R0A2  =  6  Rg2).  To  calculate  these  two  values,  which  are  used  to  describe  the 
dimensions  of  a  polymer  chain,  it  is  necessary  to  calculate  the  Kuhn  length.  The  Kuhn  length  (b) 
is  a  measure  of  the  confonnational  rigidity  and  is  the  product  of  the  contour  length  over  the 
repeating  structural  unit  (10)  and  the  characteristic  ratio  (Co,).  Here,  the  characteristic  ratio  is  a 
measure  of  short-range  interactions  along  with  the  contour  length  over  the  repeating  structural 
unit.  This  parameter  can  be  calculated  through  molecular  modeling  or  obtained  through  a 
polymer  handbook.  The  radius  of  gyration  or  the  end-to-end  distance  can  then  be  calculated 
using  Rg2  =  b  Nb  /6  or  R0A2  =  b  Nb,  respectively,  where  Nb  is  the  number  of  beads 
representing  our  polymer  chain. 

There  are  several  ways  to  calculate  the  Flory-Huggins  interaction  parameter,  x,  in  this  system. 
One  way  is  to  use  a  temperature  and  composition-dependent  interaction  parameter:  Xij  =  A 
+  B/T  +  C*  density i(r).  If  this  form  of  the  interaction  parameter  is  used,  a  minimum  of  the  A 
parameter  must  be  defined.  The  values  of  some  polymer-polymer  interactions  can  be  found  in 
polymer  handbooks.  Another  way  to  define  the  interactions  is  through  solubilities:  Xij  =  V0(5; 

-  5j)2  /RT.  These  solubility  parameters  (5)  can  be  found  in  polymer  handbooks  as  well  as  from 
molecular  dynamics  simulations.  In  both  of  these  methods,  it  is  proper  to  define  VO,  which  is  a 
reference  volume  in  units  of  cm  /mol,  and  the  temperature  in  kelvins.  If  VO  is  not  defined,  it  is 
set  to  the  default  value  of  7 1 . 1  cm  /mol.  Interactions  between  species  are  scaled  by  N,  the  index 
of  polymerization,  since  the  interaction  increases  with  the  size  of  the  polymer.  Values  of  xN 
calculated  through  these  two  methods  tend  to  be  very  large,  which  may  lead  to  numerical 
instabilities.  If  this  is  the  case,  the  interactions  should  be  capped  to  a  maximum  value.  Again, 
diffusive  dynamics  can  also  be  used  in  the  same  way  as  in  the  theoretical  case.  Two  different 
example  codes  for  an  experimental  system  composed  of  a  styrene  polymer  and  a  benzene  solvent 
are  shown  in  codes  2a  and  2b.  See  section  6  for  an  example  of  how  to  calculate  some  of  these 
experimental  values. 
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Code  2a:  Sample  Code  for  Theoretical 
species 
{ 

VO  =  70.0; 
temp  =  500; 

Styrene 

{ 

type  =  gaussian; 

R0  A  2  /  M  =  0.434; 
mass  density  =  0.969; 
solubility  parameter  =  19; 

} 

Benzene 

{ 

type  =  solvent; 
mass  density  =  0.874; 

M0  =  78; 

solubility  parameter  =  18.8; 

} 

chi 

{ 

max  chi*<N>  =  90; 
StyreneBenzene 
{ 

solubility  parameter  =  true; 

} 

} 

} 


Code  2b:  Sample  Code  for  Experimental 
species 
{ 

V0  =  70.0; 
temp  =  500; 

Styrene 

{ 

Rg  A  2  /  M  =  0.134; 
mass  density  =  0.969; 
type  =  gaussian; 

}  ' 

Benzene 

{ 

type  =  solvent; 
mass  density  =  0.874; 

M0  =  78; 

} 

chi 

{ 

max  chi*<N>  =  90; 
StyreneBenzene 
{ 

A  =  0.2; 

B  =  5.0; 

} 

} 

} 


2.3  Components 

The  Component  tab  is  used  to  define  the  architecture  of  the  system.  The  basic  structure  of  this 
group  is  as  follows: 

components! 

Ns_per_N0  =  100;  \\  resolution  along  the  polymer  chain 

firstcomponentname  { 

<component  properties> 
blocks! 

block  1 ! 

<  block  1  properties> 

} 

. . .  more  blocks  . . . 

} 

subchains! 
subchain  1  ! 
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<  subchain  1  properties> 
blocks  .... 

...  subchains  ... 

} 

. . .  more  subchains  . . . 

} 

. . .  more  components  . . .  { 

} 

} 

} 

Here,  the  blocks  and  subchains  groups  are  used  to  define  the  architecture  of  the  polymer,  and  the 
block  subgroup  within  the  blocks  group  is  used  to  define  which  polymers  are  present  when  one 
transverses  down  the  polymer  backbone.  The  subchain  group  adds  chains  to  the  end  of  the  last 
polymer  block,  where  the  subchain  subgroup  can  be  used  to  add  additional  blocks.  Again,  there 
are  theoretical  and  experimental  formats  for  this  group. 

2.3.1  Theoretical  and  Experimental  Variables 

For  the  theoretical  variables,  we  will  look  at  an  AB  three-arm  diblock  copolymer  star  polymer, 
where  each  arm  is  a  symmetric  diblock  copolymer  (see  figure  2a)  in  a  C  solvent.  The  inner  star 
is  composed  of  the  B  block.  For  the  experimental  variable,  we  look  at  a  styrene-isoprene  four- 
arm  star  polymer.  In  this  example,  one  arm  is  a  styrene-isoprene  diblock  copolymer  where  the 
isoprene  end  is  attached  to  three  styrene  arms  (see  figure  2b). 


2.3.2  Theoretical  and  Experimental  Setup 

In  the  theoretical  setup,  the  volume  fraction  (f)  for  each  species  is  specified  to  be  0.5  relative  to 
each  chain,  where  the  index  of  polymerization  (N)  is  specified  to  be  100  for  each  arm.  In  this 
gel,  the  volume  fraction  of  the  star  is  0.9  and  the  solvent  is  0.1.  In  the  experimental  setup,  the 
volume  fraction  f  is  ignored,  and  instead,  the  total  molecular  weight  of  the  polymer  is  used  along 
with  the  mole  fraction  of  each  block  relative  to  the  entire  star  polymer.  In  this  setup,  only  the 
molecular  weight  of  the  first  block  of  the  main  chain  is  read  in  (Mn  =  Total  Mn  *  mole  fraction 
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[first  block]),  and  the  molecular  weight  of  the  additional  blocks  are  calculated  internally.  In  the 
experimental  input,  it  is  not  necessary  to  calculate  N;  instead,  the  program  will  calculate  Nusing 
Mn.  In  the  experimental  setup,  it  is  also  possible  to  input  the  molecular  weight  of  each  block  and 
the  number  of  beads  that  represents  each  block  (Ns)  instead  of  the  number  of  beads  used  to 
represent  the  entire  polymer  (Ns_per_N0).  An  example  code  for  the  theoretical  and 
experimental  system  can  be  found  in  codes  3a  and  3b,  respectively. 


Code  3a:  Sample  Code  for  Theoretical 
components  { 

Ns_per_N0  =  100; 
star  { 

type  =  polymer; 
volfrac  =  0.9; 

N=  100; 
blocks  { 

block  1  { 

species  =  A; 
f  =  0.5; 

} 

block  2  { 

species  =  B; 
f=0.5; 

} 

} 

subchains  { 
subchain  1  { 

N=  100; 
blocks { 
block  1  { 
species  =  B; 
f  =  0.5; 

} 

block  2  { 
species  =  A; 
f  =  0.5; 

} 

} 

subchain  2  { 

N=  100; 
blocks { 
block  1  { 
species  =  B; 
f  =  0.5; 

} 

_ block  2  { _ 


Code  3b:  Sample  Code  for  Experimental 
components] 

Ns_per_N0  =  30; 
diblock  { 
type  =  polymer; 
volfrac  =  1.0; 
blocks] 
styrene  block] 
species  =  Styrene; 
mole  fraction  =  0.2; 

Mn  =  40000; 

} 

} 

subchains] 
styrene  subchain] 
blocks] 

A] 

species  =  Styrene; 
mole  fraction  =  0.2; 

} 

} 

} 

isoprene  subchain] 
blocks] 

A] 

species  =  Isoprene; 
mole  fraction  =  0.2; 

} 

} 

subchains] 
styrene  subchain  1  ] 
blocks ] 

A] 

species  =  Styrene; 
mole  fraction  =  0.2; 

} 

} 

I _ 
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species  =  A; 

styrene  subchain  2  { 

f  =  0.5; 

blocks { 

} 

A{ 

} 

species  =  Styrene; 

} 

mole  fraction  =  0.2; 

} 

} 

} 

} 

solvent! 

} 

type  =  solvent; 

} 

volfrac  =  0.1; 

} 

species  =  C; 

} 

} 

} 

} 

} 

3.  Defining  Simulation  Methods  and  Parameters 


Now  that  the  polymer  system  is  defined,  we  will  go  over  some  simulation  methods,  such  as 
numerical  methods,  format  and  intervals  of  output  files,  initial  conditions,  and  iteration  control. 

3.1  Initial  Field  Condition 

This  is  typically  the  first  line  in  the  input  file  and  it  specifies  the  initial  format  of  the  field.  When 
the  initial  guess  property  is  set  as  initialguess  =  random.  The  fields  are  produced  using  a 
random  number  generator.  For  debugging,  it  is  easier  to  use  a  constant  zero  initial  condition, 
which  can  be  specified  by  setting  initial  guess  =  zero.  If  field  configurations  are  appended  to  the 
end  of  the  input  file,  they  will  be  read  in  by  setting  initial  guess  =  attached,  or  if  the  initial 
condition  is  written  to  a  separate  file,  the  initial  guess  is  set  to  the  file  name.  For  example,  if  the 
fields  are  saved  in  a  file  named  fields.txt,  they  can  be  read  in  by  specifying  initial_guess  = 
fields.txt.  It  is  also  possible  to  use  spectral  seeding  to  guess  the  symmetry  of  the  system.  This 
will  be  covered  later  in  this  report. 

3.2  Iteration  Control 

This  group  specifies  the  maximum  number  of  iterations,  maxiter,  that  will  be  conducted  if  the 
minimum  error,  min  error,  is  not  reached.  If  the  minimum  error  is  reached  before  the  maximum 
iteration,  the  simulation  will  finish.  Within  this  subgroup  it  is  also  possible  to  ramp  multiple 
parameters  such  as  X,  angles,  or  simulation  box  size  length  after  a  specified  iteration  count  is 
met.  It  is  also  optional  to  specify  how  the  parameter  is  sampled,  either  linear  or  log.  If  it  is  not 
specified,  linear  sampling  is  assumed.  The  iteration  control  group  is  defined  as  follows  for  a 
maximum  of  1000  iterations  and  a  minimum  error  of  le-5  for  seven  box  sizes.  The  x  box  length 
(d_l)  is  varied  from  8  to  14  (abs_range),  by  adding  1  to  the  length  6  times  (n_step). 


8 


iteration  control  { 
max_inter  =  1000; 
min_error  =  l.e-5; 
vary  { 

1  { 

parameter  =  d_l; 
abs_range  =  8.,  14.; 
nstep  =  6; 
first_step  =  1 ; 

} 

} 

} 

3.3  Field  Relaxation 

To  update  the  fields  in  self-consistent  field  theory  through  optimization  or  diffusive  dynamics 
and  in  complex  Langevin  theory,  it  is  necessary  to  specify  desired  numerical  schemes  and  the 
time  steps. 

For  diffusive  dynamics,  the  scheme  should  be  set  to  diffusive  dynamics,  and  the  diffusive  time 
step  should  be  specified.  Here,  the  fields  are  updated  every  iteration  (n  iter),  starting  from  the 
first  iteration  (first  iter).  Here,  the  sis  const  is  the  relaxation  constant  for  the  pressure  field 
update. 

field  relaxation 

{ 

firstiter  =  0; 
niter  =  1 ; 

scheme  =  diffusive  dynamics; 
sis_const  =  5.0 
diffusive  time  step=l; 

} 

When  self-consistent  field  theory  optimization  or  complex  Langevin  is  used,  there  are  two 
methods  that  can  be  used:  the  explicit  Euler  method  or  the  semi-implicit-seidel  method. 

The  first  method  is  the  explicit  Euler  method.  When  scft  optimization  is  used,  this  group  takes 
the  following  fonn  where  euler  const  l  is  the  relaxation  constant  for  the  pressure  field  and 
euler_const_2  is  the  relaxation  constant  for  the  nonpressure  field. 

field  relaxation 

{ 

firstiter  =  0; 
niter  =  1 ; 
scheme  =  euler; 
euler_const_l  =  1; 
euler_const_2  =  0.1; 
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} 


If  complex  Langevin  is  instead  desired,  it  is  necessary  to  add  two  additional  parameters  to  the 
group:  complex  langevin  =  true  and  ginzburg  =100. 

Another  more  stable  method  to  update  the  fields  is  the  semi-implicit-seidel  method  developed  by 
Ceniceros  and  Fredrickson  ( 4 ).  Here  sis  const  is  the  field  relaxation  constant  for  the  pressure 
field,  while  sisconstminus  is  the  relaxation  constant  for  the  nonpressure  field.  It  is  also 
possible  to  allow  the  value  of  the  relaxation  constants  to  be  rescaled  through  the  simulation  by 
including  auto  lambda  =  true  and  auto  lambda  restart  =  true  (can  also  be  used  in  diffusive 
dynamics  and  in  the  Euler  method). 

field  relaxation 

{ 

firstiter  =  0; 

niter  =  1 ; 

scheme  =  sis; 

sis_const  =  5.; 

sisconstminus  =  1.; 

auto  lambda  =  true; 

auto  lambda  restart  =  true; 

} 

3.4  Status  Save 

If  there  is  a  parameter  value  that  you  would  like  calculated  every  n  iter  iterations,  such  as  the 
free  energy,  box  size,  and  angles,  the  Status  Save  group  can  be  used  to  store  the  output.  Each 
variable  that  you  would  like  to  save  must  be  on  a  new  line  followed  by  the  variable  type  and  the 
type  of  output. 

In  the  following  example,  the  code  will  save  data  to  a  file  at  every  iteration  and  will  dump  the 
iteration  number  (the  box  size  in  the  1 -direction)  and  the  Xac  into  a  file. 

status  save  { 
niter  =  1 ; 

outputs  =  restart  file->  iteration,  int,  %05d 
results->F,  double,  %24.16g 
box  size->d_l,  double, %8.3f; 
species->  chi->A_C->A, double, %4.2f 

} 

3.5  Field  Save 

When  running  your  job,  it  is  sometimes  important  to  restart  your  job.  In  order  to  do  this,  it  is 
necessary  to  save  your  fields  throughout  the  simulation.  The  Field  Save  group  can  be  used  to 
save  the  fields  at  a  specified  frequency  or  after  the  error  change  is  greater  than  a  given  threshold. 
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In  addition  to  the  restart  capability,  the  data  saved  using  this  group  can  also  be  used  to  visualize 
results  before  the  simulation  is  completed. 

In  the  following  example,  the  fields  are  saved  in  text  format  every  100  iterations  or  if  the  error 
change  is  >0.4,  but  not  after  each  iteration. 

field  save{ 

niter  =  100; 
format  =  text; 
save  recent  =  false; 
delta  error  =  0.4; 

} 

For  restarting  the  run,  the  data  should  be  saved  in  text  fonnat,  but  for  visualization,  two 
additional  formats  can  be  saved:  gnuplot  and  binary. 

3.6  Misc 

This  group  dictates  the  numerical  methods  used  to  solve  the  modified  diffusion  equation  for  the 
propagators  and  the  quadrature  used  to  calculate  the  densities.  In  addition,  the  verbosity  of  the 
simulation  is  indicated  through  this  group.  The  verbose  level  indicates  how  much  information  is 
sent  to  the  tenninal  or  to  the  standard  output,  where  a  level  of  0  will  produce  no  output  and  a 
level  of  6  will  print  out  the  value  of  every  array  at  every  iteration.  High  verbose  levels  are  used 
for  debugging  purposes. 

Two  methods  used  to  calculate  the  modified  diffusion  equation  are  the  operator  splitting  method 
(op  split)  or  the  backward  differentiation  formula  (bdf)  (5).  Two  quadratures  can  be  used, 
Gaussian  (gauss)  or  Simpson  (simpson). 

misc  { 

verbose  level  =  2; 
quadrature  =  gauss; 
diffusion  equation  { 

algorithm  =  opsplit; 

} 

} 


4.  Extensions 


4.1  Unit  Cell  Calculations 

To  calculate  an  optimal  unit  cell  of  your  morphology,  you  can  use  the  box  relaxation  group, 
which  will  resize  the  simulation  box  length  and/or  box  angles  to  an  optimal  cell  with  minimum 
stress. 
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In  the  following  example,  the  box  is  resized  every  iteration  after  the  fifth  iteration  as  long  as  the 
error  is  <l.e-3  and  the  total  stress  on  the  box  is  >l.e-4.  The  length  d_l  is  relaxed  at  a  rate 
constant  of  0.1,  where  d_2  is  constrained  to  the  same  length.  The  third  length  d_3  and  all  the 
angles  are  not  allowed  to  change. 

box  relaxation! 
first_iter  =  5; 
niter  =  1 ; 
min_error=  l.e-3; 
min_stress  =  l.e-4; 
speed  =  0.1; 
which_dims{ 
d_l  =  true; 
d_2  =  constrain(d_l); 
d_3  =  false; 
ang_l_l  =  false; 
ang_l_2  =  false; 
ang_2_3  =  false; 
d_l  bounds  =  10,15; 
d_2  bounds  =  10,15; 

} 

} 

4.2  Spectral  Seeding 

It  is  possible  to  seed  in  simple  symmetries  into  the  simulation.  One  can  use  centrosymmetric 
cosine  seeding  (-id)  or  plane  wave  basis  seeding  (id). 

spectral) 

space  group  =  -id; 
weights  =  4,  1; 

} 

The  weights  indicate  a  prefactor  for  each  harmonic.  For  this  case,  only  the  first  two  harmonics 
are  used:  4  cos(27tx/D)  +  cos(47tx/D).  If  more  weights  are  included,  more  harmonics  are  used. 

4.3  Ensemble  of  Particles 

To  implement  particles,  we  can  use  the  ensemble  group,  where  the  nanoparticles  will  be  moved  a 
total  of  specified  moves  (totalmoves).  The  particles  will  be  moved  after  the  error  tolerance 
specified  in  the  Iteration  Control  group  is  reached.  In  setting  up  the  nanoparticles,  a  specified 
number  of  position  updates  will  be  performed  between  additions  of  particles  (moves_per_step). 
The  ratio  of  the  variable  delta  over  the  variable  ginzburg  in  the  Field  Relaxation  group  defines 
the  ratio  between  the  force-biased  tenn  and  the  random  term  in  the  position  update,  where  delta 
is  the  step  width  that  specifies  the  maximum  distance  each  particle  can  move  at  a  given  step. 

The  method  variable  indicates  the  method  that  will  be  used  to  move  the  particles.  There  are 
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three  choices:  smart  Monte  Carlo  (SMC),  Langevin  dynamics  (langevin),  and  standard 
Metropolis  Monte  Carlo  (MC).  The  three  methods  are  performed  slightly  differently.  In  SMC, 
the  particles  are  displaced  based  on  the  force  and  a  random  force  based  on  the  step  width. 
Langevin  displaces  particles  using  the  same  method  as  SMC  where  all  moves  are  accepted  unless 
it  leads  to  overlap,  but  the  step  size  must  be  small  for  the  calculation  to  be  accurate  and  stable. 
MC  displaces  particles  based  on  a  random  force  calculated  from  the  step  width,  but  it  could  lead 
to  artificial  clustering  of  particles.  The  maximum  distance  a  particle  can  move  is  specified  using 
the  cutoff  variable.  The  thickness  and  overlapmargin  variables  specify  the  width  of  the 
interfacial  region  of  the  particle  cavity  function  and  the  margin  used  to  prevent  particles  from 
overlapping,  respectively. 

In  the  following  example,  the  4  particles  of  radius  75  are  moved  100  times  where  they  are 
initially  placed  at  random  with  one  position  step  between  particle  additions.  The  particles  are 
allowed  to  move  a  maximum  step  width  of  125  using  the  SMC  algorithm.  Thickness  and 
overlap  margins  of  25  are  used. 

ensembles 

{ 

total_moves=  100; 
moves_per_step  =  1 ; 
delta  =125; 
cutoff  =  1000; 
method  =  SMC; 
setupcondition  =  random; 
rotations_per_move  =  0; 
maxrotation  =  0.0; 
freq_reinsert  =  0; 
numreinsert  =  0; 
ensemble  1 
{ 

species  =  nano; 
radius  =  75; 
thickness  =  25.000; 

overlapmargin  =  25.000;  totalnumber  =  4; 
add_per_step  =  4; 
type  =  sphere; 

} 

} 


5.  Visualization 


To  visualize  this  data,  it  is  possible  to  use  a  variety  of  programs.  For  2-D  data,  gnuplot  is  the 
easiest  to  use  since  it  is  possible  to  have  the  program  output  the  data  in  this  fonnat.  For  3-D 
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visualization,  there  are  a  variety  of  other  programs  that  can  be  used  to  visualize  data,  such  as 
Matlab,  IFrIT  ( 6 ),  and  a  viewer  produced  by  Eric  Cochran  specifically  for  this  program.  All 
three  of  these  programs  will  involve  saving  data  in  the  text  format,  but  additional  scripts  or 
programs  must  be  used  to  extract  the  field  values  from  these  files  when  IFrIT  or  Matlab  are  used. 
For  these  scripts,  please  email  tanya.chantawansri@arl.army.mil. 


6.  Conversion  of  Experimental  Data  for  Mesotek+ 


We  will  consider  a  poly(styrene-b-isoprene)  in  the  solvent  tetradecane  (C-14).  The  phase 
diagram  can  be  found  in  Bang  et  al.  (7). 

The  experiments  should  provide  the  molecular  weight/total  mass  of  each  polymer  species  in  the 
block  copolymer  or  the  total  molecular  weight  and  mole/weight/mass  fractions. 

In  this  report,  we  are  given  the  mass  of  styrene  and  isoprene  in  the  diblock  copolymer: 

Diblock:  mass  of  styrene  =  15,200  g/mol  and  mass  of  isoprene  =  15,400  g/mol. 

From  these  species,  we  can  calculate  the  physical  properties  of  the  monomers  (see  table  1). 
These  properties  were  obtained  via  the  program  Synthia  (5)  of  Material  Studio,  but  they  can  also 
be  found  in  polymer  handbooks.  For  the  solvent,  the  repeat  unit  and  characteristic  ratio  are  not 
required. 


Table  1.  Properties  for  isoprene,  styrene,  and  C-14  obtained  using  Synthia. 


Species 

Monomer  Volume 
(cm3/mol) 

Monomer  Mass 
(g/mol) 

Repeat  Unit 
Length 

(A) 

Characteristic 

Ratio 

Solubility 

Parameter 

Isoprene 

76.6 

68.1 

4.4 

5.4 

16.9 

Styrene 

97.0 

104.2 

2.5 

9.9 

20.1 

C-14 

225.1 

196.3 

NA 

NA 

17.5 

Note:  NA  =  not  applicable. 


First  it  is  necessary  to  calculate  the  volume  fraction  of  isoprene  and  styrene.  This  can  be  done 
using  the  density  of  each  component:  styrene  density  =  1.047  g/crn  and  isoprene  density 
=  0.913  g/cm3. 

The  volume  of  styrene  is  calculated  by  dividing  the  mass  of  styrene  in  the  diblock  by  the  density 
of  styrene:  (15,200  g/mol)/(1.047  g/cm3)  =  14517.67  cnrVmol. 

The  volume  of  isoprene  is  calculated  by  dividing  the  mass  of  isoprene  in  the  diblock  by  the 
density  of  isoprene:  (15,400  g/mol)/(0.913  g/cnr’)  =  16,867.5  cm3/mol. 
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The  volume  fraction  of  styrene  (Fs)  and  isoprene  (Fi)  in  the  diblock  can  then  be  calculated  by 
dividing  the  volume  of  the  species  of  interest  by  the  total  volume: 

Fs  =  14517.67  cm3/mol  /  (14517.67  cm3/mol  +  16867.5  cm3/mol)  =  0.46. 

Fi=  16867.5  cm3/mol/  (14517.67  cm3/mol  +  16867.5  cm3/mol)  =  0.54. 

Using  these  values,  we  can  calculate  the  volume  and  molecular  weight  of  a  bead.  Let  us  define 
our  polymer  as  having  60  total  beads,  which  will  give  us  28  styrene  beads  and  32  isoprene  beads. 

Relevant  equations  include  the  total  molecular  weight  MN  =  niMi  +  nsMs,  the  volume  fraction  of 
each  species  Fi  =  n(Vi  /(niVi  +  nsVs),  and  Fs  =  nsVs  /(niVi  +  nsVs),  where  MN  is  the  total  mass  of 
a  polymer  chain,  Mi  is  the  mass  of  an  isoprene  monomer,  Ms  is  the  mass  of  a  styrene  monomer, 

Fi  is  the  volume  fraction  of  an  isoprene  monomer,  Fs  is  the  volume  fraction  of  a  styrene 
monomer,  Vi  is  the  volume  occupied  by  a  monomer  of  isoprene,  and  Vs  is  the  volume  occupied 
by  a  monomer  of  styrene. 

Combining  these  three  equations,  we  define  a  new  parameter  C:  C  =  MN  /  (Mi  Fi  Vs  +  Ms  Fs 
Vi).  Using  the  values  calculated  through  Synthia  (5),  we  get  the  following  value  for  C: 

C  =  (15,200  g/mol  +  15,400  g/mol)  /  (68.1g/mol  x  0.54  x  97.0  cm3/mol  +  104.2  g/mol  x  0.46 
x  76.6  cmVmol)  =  4.23  mol/cm3. 

We  can  then  calculate  the  number  of  isoprene  and  styrene  monomers  denoted  by  ns  and  ni, 
respectively. 

ni  =  C  Fi  Vs=  4.23  mol/cm3  x  0.54  x  97. 0  cm3/mol  =  221  monomers, 
ns  =  C  Fs  Vi  =  4.23  mol/cm3  x  0.46  x  76.6  cm3/mol  =  149  monomers. 

We  then  need  to  calculate  how  many  monomers  are  in  a  bead  for  each  component.  For  styrene, 
we  have  149  monomers  in  28  beads  =  5.32.  For  isoprene,  we  have  221  monomers  in  32  beads 
=  6.91. 

To  calculate  the  volume  of  each  bead,  we  multiply  the  volume  for  a  monomer  by  the  following 
ratio: 

Bead  volume  for  styrene  =  5.32  x  97. 0  cm  /mol  =  516.04  cm  /mol. 

Bead  volume  for  isoprene  =  6.91  x  76.6  cm  /mol  =  529.31  cm  /mol. 

To  get  the  mass  of  each  bead,  we  multiply  the  mass  of  a  monomer  by  the  following  ratio: 

Bead  mass  of  styrene  =  5.32  x  104.2  g/mol  =  554.34  g/mol. 

Bead  mass  of  isoprene  =  6.91  x  68.1  g/mol  =  470.57  g/mol. 

To  get  the  Kuhn  length,  we  multiply  the  repeat  unit  length  and  characteristic  ratio: 

Kuhn  length  for  styrene  =  2.5  A  x  9.9  =  24.8  A. 

Kuhn  length  for  isoprene  =  4.4  A  x  5.4  =  23.8  A. 
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The  Kuhn  length  is  related  to  the  radius  of  gyration.  To  calculate  the  parameter  RgA2/M,  we 
multiply  the  square  of  the  Kuhn  length  and  divide  by  the  product  of  the  bead  mass  and  the 
number  6: 

RgA2/M  for  styrene  =  24.8  A  x  24.8  A  /  (6  x  554.3  g/mol)  =  0.185  A2g/mol. 

RgA2/M  for  isoprene  =  23.8  A  x  23.8  A  /  (6  x  470.6  g/mol)  =  0.201  A2g/mol. 

Since  we  are  representing  the  solvent  as  one  bead,  we  can  use  the  parameters  obtained  by 
Synthia  (5)  for  our  calculations.  Table  2  summarizes  the  parameters  that  should  be  inputted  into 
Mesotek+. 

Table  2.  Calculate  input  parameters  for  Mesotek+. 


Species 

Bead  Volume 
(cm3/mol) 

Bead  Mass 
(g/mol) 

Solubility 

RgA2/M 

(A2g/mol) 

Styrene 

516.0 

554.3 

20.1 

0.185 

Isoprene 

529.3 

470.6 

16.9 

0.201 

C14 

225.1 

196.3 

17.5 

NA 

Note:  NA  =  not  applicable. 


The  following  pages  present  a  script  for  this  experimental  system.  Here,  the  molecular  weight  of 
each  block  is  off  from  what  is  given  experimentally  due  to  roundoff  errors  when  the  bead  mass 
and  volume  were  calculated,  but  this  difference  should  not  affect  the  results.  Here,  a  2-D 
simulation  of  a  40/60  poly(styrene-b-isoprene)/C14  mixture  at  373  K  is  considered,  which 
should  correspond  to  the  cylindrical  phase.  Results  are  shown  in  figure  3. 


Figure  3.  Results  from  40/60  volume  styrene-b-isoprene  +  tetradecane  using  gnuplot:  A)  styrene  volume  fraction, 
B)  isoprene  volume  fraction,  and  C)  tetradecane  volume  fraction.  The  color  bar  to  the  right  of  each  plot 
indicates  how  the  coloring  corresponds  to  the  volume  fraction. 


initial  guess  =  random; 

components 

{ 

polymer 

{ 

type  =  polymer; 
volfrac  =  0.4; 
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blocks 

{ 

Bead  Styrene 

{ 

species  =  Styrene; 

Mn  =  15520.4;  //Ns  *  Bead  Mass  ~  152,000 
Ns  =  28; 

} 

Bead  Isoprene 

{ 

species  =  Isoprene; 

Mn  =  1 5059.2;  //Ns  *  Bead  Mass  ~  1 54,000 

Ns  =  32; 

} 

} 

} 

cl4 

{ 

type  =  solvent; 
volfrac  =  0.6; 
species  =  C14; 

} 

} 

species 

{ 

temp  =  373; 

C14 

{ 

type  =  solvent; 
solubility  parameter  =  17.5; 

M0  =  196.3; 

Rg2/M  =  0.217219; 
molar  volume  =  225.1 ; 

} 

Styrene 

{ 

type  =  gaussian; 
solubility  parameter  =  20.1; 

M0  =  554.3; 

Rg2/M  =  0.18493; 
molar  volume  =  516; 

} 

Isoprene 

{ 

type  =  gaussian; 
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solubility  parameter  =  16.9; 
MO  =  470.6; 

Rg2/M  =  0.201; 
molar  volume  =  529.3; 

} 

chi 

{ 

max  chi*<N>  =  90; 
C14_Styrene 
{ 

solubility  parameter  =  true; 

} 

C14_Isoprene 

{ 

solubility  parameter  =  true; 

} 

IsopreneStyrene 

{ 

solubility  parameter  =  true; 

} 

} 

} 

iteration  control 

{ 

max_iter=  10000; 
min_error_plus  =  0.0001; 
minerrorminus  =  0.0001; 

} 

field  relaxation 

{ 

first  iter  =  1 ; 
niter  =  1 ; 
scheme  =  sis; 
plus  scheme  =  sis; 
minus  scheme  =  sis; 
sis_const  =  5; 
sis_const_minus  =  0.2; 
auto  lambda  =  true; 

} 

box  size 

{ 

rank  =  2; 
n_l =  32; 
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n_2  =  32; 
d_l =  320; 
d_2  =  320; 
ang_l_2  =  90; 

} 

field  save 

{ 

niter  =100; 
format  =  gnuplot; 
delta  =  50; 
save  recent  =100; 

} 


misc 

{ 

quadrature  =  simpson; 

diffusion  equation 

{ 

algorithm  =  split; 
extra  sampling  =  1 ; 

} 

} 
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